
### Preliminaries
source("psrm_scriptX.R")

library(Hmisc)
library(xtable)
library(lme4)
library(dynsim)
library(mvtnorm)

#set.seed(12345)

## 2.2 Descriptive Statistics
macro$yearf = as.factor(macro$year)
macro$countryf = as.factor(macro$country_name)

# Table A1
foe = macro[ , c("market" , "direct","top_incrate_n", "topitaxrate2", "l1.tradeallyadv" ,"l1rural.ineq", "l1.e2","votetax2","l1.rgdppc", "l1.trade")]
stargazer(foe, digits = 2, type = "text")

### Run models to generate coverage counts excluding any missing variables
market.empty.d = lm(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq 
                    + l1.unisuffrage + votetax2 + l1.rgdppc + l1.trade 
                    + countryf + yearf, data = macro)
ii2.marketx = lm(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq 
                 + l1.unisuffrage + votetax2 + l1.rgdppc + l1.trade 
                 + countryf + year, data = macro)

c.obs = rep(NA, length(unique(market.empty.d$model$countryf)))
minyear = rep(NA, length(c.obs))
maxyear = rep(NA, length(c.obs))
for(i in 1:length(unique(market.empty.d$model$countryf))){
  c.obs[i] = length(which(market.empty.d$model$countryf == unique(market.empty.d$model$countryf)[i]))
  minyear[i] = min(ii2.marketx$model$year[which(market.empty.d$model$countryf == unique(market.empty.d$model$countryf)[i])])
  maxyear[i] = max(ii2.marketx$model$year[which(market.empty.d$model$countryf == unique(market.empty.d$model$countryf)[i])])
}

c.ob.frame = as.data.frame(unique(market.empty.d$model$countryf))
names(c.ob.frame) = "Country"
c.ob.frame$Country = as.character(levels(c.ob.frame$Country))[c.ob.frame$Country]
c.ob.frame$Country = capitalize(c.ob.frame$Country)
c.ob.frame$Country[which(c.ob.frame$Country == "Uk")] = "United Kingdom"
c.ob.frame$Observations = c.obs
c.ob.frame$First = minyear
c.ob.frame$Last = maxyear

c.ob.frame


## 2.3 Robustness: State Capacity

market.preferred.sc1 = lm(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                          + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.soldiers
                          + as.factor(country_name) + as.factor(year), data = macro)

direct.preferred.sc1 = lm(direct ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                          + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.soldiers
                          + as.factor(country_name) + as.factor(year), data = macro)

ytax.preferred.sc1 = lm(top_incrate_n ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                        + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.soldiers
                        + as.factor(country_name) + as.factor(year), data = macro)

htax.preferred.sc1 = lm(topitaxrate2 ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                        + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.soldiers
                        + as.factor(country_name) + as.factor(year), data = macro)

market.preferred.sc2a = lm(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                           + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.edu 
                           + as.factor(country_name) + as.factor(year), data = macro)

direct.preferred.sc2a = lm(direct ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                           + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.edu 
                           + as.factor(country_name) + as.factor(year), data = macro)

ytax.preferred.sc2a = lm(top_incrate_n ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                         + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.edu 
                         + as.factor(country_name) + as.factor(year), data = macro)

htax.preferred.sc2a = lm(topitaxrate2 ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                         + l1.e2 + l1.trade + l1.rgdppc + votetax2 +  l1.edu 
                         + as.factor(country_name) + as.factor(year), data = macro)

market.preferred.sc2b = lm(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                           + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.edu + l1.jai
                           + as.factor(country_name) + as.factor(year), data = macro)

direct.preferred.sc2b = lm(direct ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                           + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.edu + l1.jai
                           + as.factor(country_name) + as.factor(year), data = macro)

ytax.preferred.sc2b = lm(top_incrate_n ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                         + l1.e2 + l1.trade + l1.rgdppc + votetax2 + l1.edu + l1.jai
                         + as.factor(country_name) + as.factor(year), data = macro)

htax.preferred.sc2b = lm(topitaxrate2 ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                         + l1.e2 + l1.trade + l1.rgdppc + votetax2 +  l1.edu + l1.jai
                         + as.factor(country_name) + as.factor(year), data = macro)
### Table A2
stargazer(market.preferred.sc1, direct.preferred.sc1, ytax.preferred.sc1, htax.preferred.sc1, omit = c("year", "country", "Constant")
          #, type = "text"
          , order = c("tradeallyadv", "rural.ineq", "lta.ineq", "rgdppc", "votetax", "e2", "trade", "taxgdp")
          , dep.var.labels = c("Indirect Share", "Direct Share", "Top Income", "Top Inheritance")
          , covariate.labels = c("Labour trade advantage (LTA)", "Inequality", "LTA : inequality", "GDP per capita", "Vote-tax link","Economic franchise",  "Trade",  "Military capacity")
          , style = "ajps", omit.stat = c("ser", "f"), add.lines = list(c("Country fixed effects", rep("Y", 4)), c("Year fixed effects", rep("Y", 4))), digits = 2
          , table.placement = "p", title = "The determinants of tax progressivity, 1870-1913: robustness to state capacity control (1)", label = "tab:robCapacity"
          , type = "text"
)


### Table A3

stargazer(market.preferred.sc2a, market.preferred.sc2b, direct.preferred.sc2a, direct.preferred.sc2b,  ytax.preferred.sc2a, ytax.preferred.sc2b, htax.preferred.sc2a, htax.preferred.sc2b, omit = c("year", "country", "Constant")
          #, type = "text"
          , order = c("tradeallyadv", "rural.ineq", "lta.ineq", "rgdppc", "votetax", "e2", "trade", "edu", "jai")
          , dep.var.labels = c("Indirect Share", "Direct Share", "Top Income", "Top Inheritance")
          , covariate.labels = c("Labour trade advantage", "Inequality", "LTA : inequality", "GDP per capita", "Vote-tax link","Economic franchise",  "Trade",  "Education", "Education : franchise")
          , style = "ajps", omit.stat = c("ser", "f"), add.lines = list(c("Country fixed effects", rep("Y", 8)), c("Year fixed effects", rep("Y", 8))), digits = 2
          , table.placement = "p", title = "The determinants of tax progressivity, 1870-1913: robustness to state capacity control (2)", label = "tab:robCapacity2"
          , type = "text"
)

### Figure A2
###### a figure for the inheritance tax results controlling for education and education-suffrage interaction (supplementary materials):


beta.hat <- coef(htax.preferred.sc2b)
# pull out the covariance matrix
cov <- vcov(htax.preferred.sc2b)
# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(htax.preferred.sc2b$model$l1rural.ineq), max(htax.preferred.sc2b$model$l1rural.ineq), length.out = 1000)
# calculate the instantaneous effect of x as z varies
dy.dx <- beta.hat["l1.tradeallyadv"] + beta.hat["l1.lta.ineq"]*z0
# calculate the standard error of each estimated effect
se.dy.dx <- sqrt(cov["l1.tradeallyadv", "l1.tradeallyadv"] + z0^2*cov["l1.lta.ineq", "l1.lta.ineq"] + 2*z0*cov["l1.tradeallyadv", "l1.lta.ineq"])
# calculate upper and lower bounds of a 95% CI
upr <- dy.dx + 1.96*se.dy.dx
lwr <- dy.dx - 1.96*se.dy.dx

# plot the ME using compactr

#par(mfrow = c(1,1))
par(mar = c(4,5,2,2)+0.4)
plot(density(htax.preferred.sc2b$model$l1rural.ineq), type = "l",  xlim = mm(htax.preferred.sc2b$model$l1rural.ineq), yaxt = "n", xaxt = "n", ylab = "", xlab = "", col = "transparent", main = "")
dens = density(htax.preferred.sc2b$model$l1rural.ineq)

mnri = min(htax.preferred.sc2b$model$l1rural.ineq, na.rm = T)
mxri = max(htax.preferred.sc2b$model$l1rural.ineq, na.rm = T)
x1 = min(which(dens$x > mnri))  
x2 = max(which(dens$x < mxri))  
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray85", border = "gray85"))

par(new = T)
plot(z0, dy.dx, lwd = 3, type = "l", xlim = mm(htax.preferred.sc2b$model$l1rural.ineq), ylim = mm(c(upr, lwr)),
     xlab = "Inequality",
     ylab = expression(frac(partialdiff*"top inheritance tax rate", partialdiff*"LTA")))
abline(h = 0, col = "grey60", lty = 2)
#lines(z0, dy.dx, lwd = 3)
lines(z0, lwr, lty = 3)
lines(z0, upr, lty = 3)

dev.print(device = pdf, file = "figureA2.pdf", width = 5, height = 5)

## 2.4 Robustness: Alternative Specifications

# recenter variables for use in HLMs:
macro$l1.rgdppcTenKs = macro$l1.rgdppc/10000
macro$l1.edupct = macro$l1.edu/100
macro$l1.jaipct = macro$l1.jai/100
macro$l1.market = macro$l1.market/100
macro$l1.direct = macro$l1.direct/100

macro$countrych = as.character(levels(macro$countryf))[macro$countryf]
macro$countrych[which(macro$countrych == "united kingdom")] = "unitedkingdom"
macro$countryf = as.factor(macro$countrych)

### Lagged dependent variable models
market.preferred.asldv = lm(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                            + l1.e2 + l1.trade + l1.rgdppcTenKs + votetax2 
                            + l1.market	+ countryf + as.factor(year), data = macro)

direct.preferred.asldv = lm(direct ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                            + l1.e2 + l1.trade + l1.rgdppcTenKs + votetax2 
                            + l1.direct + countryf + as.factor(year), data = macro)

ytax.preferred.asldv = lm(top_incrate_n ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                          + l1.e2 + l1.trade + l1.rgdppcTenKs + votetax2 
                          + l1.topincome + countryf + as.factor(year), data = macro)

htax.preferred.asldv = lm(topitaxrate2 ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                          + l1.e2 + l1.trade + l1.rgdppcTenKs + votetax2 
                          + l1.topinher + countryf + as.factor(year), data = macro)

### Random effects/HL models
market.preferred.asre = lmer(market ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                             + l1.e2 + l1.trade + l1.rgdppcTenKs + votetax2 
                             + as.factor(year) + (1 | country_name), data = macro)

direct.preferred.asre = lmer(direct ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                             + l1.e2 + l1.trade + l1.rgdppcTenKs + votetax2 
                             + as.factor(year) + (1 | country_name), data = macro)

ytax.preferred.asre = lmer(top_incrate_n ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                           + l1.e2 + l1.trade + l1.rgdppcTenKs + votetax2 
                           + as.factor(year) + (1 | country_name), data = macro)

htax.preferred.asre = lmer(topitaxrate2 ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                           + l1.e2 + l1.trade + l1.rgdppcTenKs + votetax2 
                           + as.factor(year) + (1 | country_name), data = macro)

### Table A4: LDV Models
stargazer(market.preferred.asldv,  direct.preferred.asldv,   ytax.preferred.asldv, htax.preferred.asldv, omit = c("year", "country", "Constant")
          , type = "text"
          , order = c("tradeallyadv", "rural.ineq", "lta.ineq", "rgdppc", "votetax", "e2", "trade", "edu", "jai")
          , dep.var.labels = c("Indirect Share", "Direct Share", "Top Income", "Top Inheritance")
          , covariate.labels = c("Labour trade advantage", "Inequality", "LTA : inequality", "GDP per capita (10k)", "Vote-tax link","Economic franchise",  "Trade",  "Domestic Indirect share$_{t-1}$", "Direct share$_{t-1}$", "Top income rate$_{t-1}$", "Top inheritance rate$_{t-1}$")
          , style = "ajps", omit.stat = c("ser", "f")
          , add.lines = list(c("Country fixed effects", rep("Y", 4)), c("Year fixed effects", rep("Y", 4))), digits = 2 
          , table.placement = "p", title = "The determinants of tax progressivity, 1870-1913: robustness to alternative specifications: lagged dependent variable models", label = "tab:robAltSpecsLDV"
    )

### Table A5: RE Models
stargazer(market.preferred.asre,  direct.preferred.asre,   ytax.preferred.asre, htax.preferred.asre, omit = c("year", "country", "Constant")
          , type = "text"
          , order = c("tradeallyadv", "rural.ineq", "lta.ineq", "rgdppc", "votetax", "e2", "trade", "edu", "jai")
          , dep.var.labels = c("Indirect Share", "Direct Share", "Top Income", "Top Inheritance")
          , covariate.labels = c("Labour trade advantage", "Inequality", "LTA : inequality", "GDP per capita (10k)", "Vote-tax link","Economic franchise",  "Trade")
          , style = "ajps", omit.stat = c("ser", "f")
          , add.lines = list(c("Country random effects", rep("Y", 4)), c("Year fixed effects", rep("Y", 4))), digits = 2 
          , table.placement = "p", title = "The determinants of tax progressivity, 1870-1913: robustness to alternative specifications: random effects models", label = "tab:robAltSpecsRE"
)

### Figures from LDV models -- figure A3

beta.hat.market <- coef(market.preferred.asldv)
# pull out the covariance matrix
cov.market <- vcov(market.preferred.asldv)
# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(market.preferred.asldv$model$l1rural.ineq), max(market.preferred.asldv$model$l1rural.ineq), length.out = 1000)
# calculate the instantaneous effect of x as z varies
dy.dx.market <- beta.hat.market["l1.tradeallyadv"] + beta.hat.market["l1.lta.ineq"]*z0
# calculate the standard error of each estimated effect
se.dy.dx.market <- sqrt(cov.market["l1.tradeallyadv", "l1.tradeallyadv"] + z0^2*cov.market["l1.lta.ineq", "l1.lta.ineq"] + 2*z0*cov.market["l1.tradeallyadv", "l1.lta.ineq"])
# calculate upper and lower bounds of a 95% CI
upr.market <- dy.dx.market + 1.96*se.dy.dx.market
lwr.market <- dy.dx.market - 1.96*se.dy.dx.market

# plot the ME using compactr

#par(mfrow = c(1,1))
par(mar = c(4,5,2,2)+0.4, mfrow = c(2, 2))

plot(density(market.preferred.asldv$model$l1rural.ineq), type = "l",  xlim = mm(market.preferred.asldv$model$l1rural.ineq), yaxt = "n", xaxt = "n", ylab = "", xlab = "", col = "transparent", main = "")
dens = density(market.preferred.asldv$model$l1rural.ineq)

mnri = min(market.preferred.asldv$model$l1rural.ineq, na.rm = T)
mxri = max(market.preferred.asldv$model$l1rural.ineq, na.rm = T)
x1 = min(which(dens$x > mnri))  
x2 = max(which(dens$x < mxri))  
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray85", border = "gray85"))

par(new = T)
plot(z0, dy.dx.market, lwd = 3, type = "l", xlim = mm(market.preferred.asldv$model$l1rural.ineq), ylim = mm(c(upr.market, lwr.market)),
     xlab = "Inequality",
     ylab = expression(frac(partialdiff*"indirect tax share", partialdiff*"LTA")),
     main = "Indirect tax share")
abline(h = 0, col = "grey60", lty = 2)
#lines(z0, dy.dx, lwd = 3)
lines(z0, lwr.market, lty = 3)
lines(z0, upr.market, lty = 3)

### direct
beta.hat.direct <- coef(direct.preferred.asldv)
# pull out the covariance matrix
cov.direct <- vcov(direct.preferred.asldv)
# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(direct.preferred.asldv$model$l1rural.ineq), max(direct.preferred.asldv$model$l1rural.ineq), length.out = 1000)
# calculate the instantaneous effect of x as z varies
dy.dx.direct <- beta.hat.direct["l1.tradeallyadv"] + beta.hat.direct["l1.lta.ineq"]*z0
# calculate the standard error of each estimated effect
se.dy.dx.direct <- sqrt(cov.direct["l1.tradeallyadv", "l1.tradeallyadv"] + z0^2*cov.direct["l1.lta.ineq", "l1.lta.ineq"] + 2*z0*cov.direct["l1.tradeallyadv", "l1.lta.ineq"])
# calculate upper and lower bounds of a 95% CI
upr.direct <- dy.dx.direct + 1.96*se.dy.dx.direct
lwr.direct <- dy.dx.direct - 1.96*se.dy.dx.direct


plot(density(direct.preferred.asldv$model$l1rural.ineq), type = "l",  xlim = mm(direct.preferred.asldv$model$l1rural.ineq), yaxt = "n", xaxt = "n", ylab = "", xlab = "", col = "transparent", main = "")
dens = density(direct.preferred.asldv$model$l1rural.ineq)

mnri = min(direct.preferred.asldv$model$l1rural.ineq, na.rm = T)
mxri = max(direct.preferred.asldv$model$l1rural.ineq, na.rm = T)
x1 = min(which(dens$x > mnri))  
x2 = max(which(dens$x < mxri))  
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray85", border = "gray85"))

par(new = T)
plot(z0, dy.dx.direct, lwd = 3, type = "l", xlim = mm(direct.preferred.asldv$model$l1rural.ineq), ylim = mm(c(upr.direct, lwr.direct)),
     xlab = "Inequality",
     ylab = expression(frac(partialdiff*"direct tax share", partialdiff*"LTA")),
     main = "Direct tax share")
abline(h = 0, col = "grey60", lty = 2)
#lines(z0, dy.dx, lwd = 3)
lines(z0, lwr.direct, lty = 3)
lines(z0, upr.direct, lty = 3)

### top income
beta.hat.ytax <- coef(ytax.preferred.asldv)
# pull out the covariance matrix
cov.ytax <- vcov(ytax.preferred.asldv)
# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(ytax.preferred.asldv$model$l1rural.ineq), max(ytax.preferred.asldv$model$l1rural.ineq), length.out = 1000)
# calculate the instantaneous effect of x as z varies
dy.dx.ytax <- beta.hat.ytax["l1.tradeallyadv"] + beta.hat.ytax["l1.lta.ineq"]*z0
# calculate the standard error of each estimated effect
se.dy.dx.ytax <- sqrt(cov.ytax["l1.tradeallyadv", "l1.tradeallyadv"] + z0^2*cov.ytax["l1.lta.ineq", "l1.lta.ineq"] + 2*z0*cov.ytax["l1.tradeallyadv", "l1.lta.ineq"])
# calculate upper and lower bounds of a 95% CI
upr.ytax <- dy.dx.ytax + 1.96*se.dy.dx.ytax
lwr.ytax <- dy.dx.ytax - 1.96*se.dy.dx.ytax


plot(density(ytax.preferred.asldv$model$l1rural.ineq), type = "l",  xlim = mm(ytax.preferred.asldv$model$l1rural.ineq), yaxt = "n", xaxt = "n", ylab = "", xlab = "", col = "transparent", main = "")
dens = density(ytax.preferred.asldv$model$l1rural.ineq)

mnri = min(ytax.preferred.asldv$model$l1rural.ineq, na.rm = T)
mxri = max(ytax.preferred.asldv$model$l1rural.ineq, na.rm = T)
x1 = min(which(dens$x > mnri))  
x2 = max(which(dens$x < mxri))  
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray85", border = "gray85"))

par(new = T)
plot(z0, dy.dx.ytax, lwd = 3, type = "l", xlim = mm(ytax.preferred.asldv$model$l1rural.ineq), ylim = mm(c(upr.ytax, lwr.ytax)),
     xlab = "Inequality",
     ylab = expression(frac(partialdiff*"top income tax rate", partialdiff*"LTA")), 
     main = "Top income tax rate")
abline(h = 0, col = "grey60", lty = 2)
#lines(z0, dy.dx, lwd = 3)
lines(z0, lwr.ytax, lty = 3)
lines(z0, upr.ytax, lty = 3)

## top inheritance
beta.hat.htax <- coef(htax.preferred.asldv)
# pull out the covariance matrix
cov.htax <- vcov(htax.preferred.asldv)
# a set of values of z to compute the (instantaneous)
# effect of x
z0 <- seq(min(htax.preferred.asldv$model$l1rural.ineq), max(htax.preferred.asldv$model$l1rural.ineq), length.out = 1000)
# calculate the instantaneous effect of x as z varies
dy.dx.htax <- beta.hat.htax["l1.tradeallyadv"] + beta.hat.htax["l1.lta.ineq"]*z0
# calculate the standard error of each estimated effect
se.dy.dx.htax <- sqrt(cov.htax["l1.tradeallyadv", "l1.tradeallyadv"] + z0^2*cov.htax["l1.lta.ineq", "l1.lta.ineq"] + 2*z0*cov.htax["l1.tradeallyadv", "l1.lta.ineq"])
# calculate upper and lower bounds of a 95% CI
upr.htax <- dy.dx.htax + 1.96*se.dy.dx.htax
lwr.htax <- dy.dx.htax - 1.96*se.dy.dx.htax


plot(density(htax.preferred.asldv$model$l1rural.ineq), type = "l",  xlim = mm(htax.preferred.asldv$model$l1rural.ineq), yaxt = "n", xaxt = "n", ylab = "", xlab = "", col = "transparent", main = "")
dens = density(htax.preferred.asldv$model$l1rural.ineq)

mnri = min(htax.preferred.asldv$model$l1rural.ineq, na.rm = T)
mxri = max(htax.preferred.asldv$model$l1rural.ineq, na.rm = T)
x1 = min(which(dens$x > mnri))  
x2 = max(which(dens$x < mxri))  
with(dens, polygon(x=c(x[c(x1,x1:x2,x2)]), y= c(0, y[x1:x2], 0), col="gray85", border = "gray85"))

par(new = T)
plot(z0, dy.dx.htax, lwd = 3, type = "l", xlim = mm(htax.preferred.asldv$model$l1rural.ineq), ylim = mm(c(upr.htax, lwr.htax)),
     xlab = "Inequality",
     ylab = expression(frac(partialdiff*"top income tax rate", partialdiff*"LTA")), 
     main = "Top inheritance tax rate")
abline(h = 0, col = "grey60", lty = 2)
#lines(z0, dy.dx, lwd = 3)
lines(z0, lwr.htax, lty = 3)
lines(z0, upr.htax, lty = 3)

dev.print(device = pdf, file = "figureA3.pdf", height = 7, width = 7)
dev.off()

### Long-run effects in LDV models -- figure A4
## Figure A4

## four scenarios: unequalLTA, equalLTA, unequalLTD, equalLTD
## four tax outcomes: market tax revenue share, direct tax rev share, top income rates, top inheritance rates
unequalLTAScenMkt = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T),  
                               l1rural.ineq  = quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                       #       l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T)*quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                           l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                               l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                               l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                               votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                               l1.market	= quantile(macro$l1.market, 0.5, na.rm = T)
)

equalLTAScenMkt = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T),  
                             l1rural.ineq  = quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                       #    l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T)*quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                          l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                             l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                             l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                             votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                             l1.market	= quantile(macro$l1.market, 0.5, na.rm = T)
)

unequalLTDScenMkt = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T),  
                               l1rural.ineq  = quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                      #  l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T)*quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                       l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                               l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                               l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                               votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                               l1.market	= quantile(macro$l1.market, 0.5, na.rm = T)
)												

equalLTDScenMkt = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T),  
                             l1rural.ineq  = quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                      #     l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T)*quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                             l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                             l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                             l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                             votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                             l1.market	= quantile(macro$l1.market, 0.5, na.rm = T)
)						
scenMkt4corners = list(equalLTAScenMkt, unequalLTAScenMkt, equalLTDScenMkt, unequalLTDScenMkt)

Sim4 <- dynsim(obj = market.preferred.asldv, ldv = 'l1.market',
               scen = scenMkt4corners, n = 20)     

dynsimGG(Sim4, ylab = "Domestic indirect tax share", color = 'Set1',leg.labels = c("Equal, high LTA", "Unequal, high LTA", "Equal, low LTA", "Unequal, low LTA"))

dev.print(device = pdf, file = "figureA4a.pdf", width = 7.5, height = 4)    
dev.off()

### Direct tax shares:
unequalLTAScenDir = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T),  
                               l1rural.ineq  = quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                          #     l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T)*quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                               l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                               l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                               l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                               votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                               l1.direct	= quantile(macro$l1.direct, 0.5, na.rm = T),
                               countryfgermany = 1
)

equalLTAScenDir = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T),  
                             l1rural.ineq  = quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                        #      l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T)*quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                             l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                             l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                             l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                             votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                             l1.direct	= quantile(macro$l1.direct, 0.5, na.rm = T),
                             countryfgermany = 1
)

unequalLTDScenDir = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T),  
                               l1rural.ineq  = quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                         #       l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T)*quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                               l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                               l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                               l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                               votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                               l1.direct	= quantile(macro$l1.direct, 0.5, na.rm = T),
                               countryfgermany = 1
)												

equalLTDScenDir = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T),  
                             l1rural.ineq  = quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                   #   l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T)*quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                     l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                             l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                             l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                             votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                             l1.direct	= quantile(macro$l1.direct, 0.5, na.rm = T),
                             countryfgermany = 1
)		

scenDir4corners = list(equalLTAScenDir, unequalLTAScenDir, equalLTDScenDir, unequalLTDScenDir)
test = list(unequalLTAScenDir)

Sim4dir <- dynsim(obj = direct.preferred.asldv, ldv = 'l1.direct',
                  scen = scenDir4corners, n = 20)                
dynsimGG(Sim4dir, ylab = "Direct tax share", color = 'Set1', leg.labels = c("Equal, high LTA", "Unequal, high LTA", "Equal, low LTA", "Unequal, low LTA"))

dev.print(device = pdf, file = "figureA4b.pdf", width = 7.5, height = 4)   
dev.off()

unequalLTAScenYtax = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T),  
                                l1rural.ineq  = quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                      #           l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T)*quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                                l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                                l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                                l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                                votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                                l1.topincome	= quantile(macro$l1.topincome, 0.5, na.rm = T),
                                countryfgermany = 1
                                
)

equalLTAScenYtax = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T),  
                              l1rural.ineq  = quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                      #         l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T)*quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                              l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                              l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                              l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                              votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                              l1.topincome	= quantile(macro$l1.topincome, 0.5, na.rm = T),
                              countryfgermany = 1
                              
)

unequalLTDScenYtax = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T),  
                                l1rural.ineq  = quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                       #          l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T)*quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                                l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                                l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                                l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                                votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                                l1.topincome	= quantile(macro$l1.topincome, 0.5, na.rm = T),
                                countryfgermany = 1
)												

equalLTDScenYtax = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T),  
                              l1rural.ineq  = quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                     #     l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T)*quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                              l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                              l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                              l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                              votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                              l1.topincome	= quantile(macro$l1.topincome, 0.5, na.rm = T),
                              countryfgermany = 1
)		

scenYtax4corners = list(equalLTAScenYtax, unequalLTAScenYtax, equalLTDScenYtax, unequalLTDScenYtax)
test = list(unequalLTAScenYtax)		

Sim4ytax <- dynsim(obj = ytax.preferred.asldv, ldv = 'l1.topincome',
                   scen = scenYtax4corners, n = 20)                
dynsimGG(Sim4ytax, ylab = "Top income tax rate", color = 'Set1', leg.labels = c("Equal, high LTA", "Unequal, high LTA", "Equal, low LTA", "Unequal, low LTA"))

dev.print(device = pdf, file = "figureA4c.pdf", width = 7.5, height = 4)   
dev.off()

### Top Inheritance Tax Rates

unequalLTAScenHtax = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T),  
                                l1rural.ineq  = quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                   #   l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T)*quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                     l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                                l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                                l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                                votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                                l1.topinher	= quantile(macro$l1.topinher, 0.5, na.rm = T),
                                countryfnetherlands = 1
                                
)

equalLTAScenHtax = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T),  
                              l1rural.ineq  = quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                  #    l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.9, na.rm = T)*quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                     l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                              l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                              l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                              votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                              l1.topinher	= quantile(macro$l1.topinher, 0.5, na.rm = T),
                              countryfnetherlands = 1
)

unequalLTDScenHtax = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T),  
                                l1rural.ineq  = quantile(macro$l1rural.ineq, 0.9, na.rm = T),
               #     l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T)*quantile(macro$l1rural.ineq, 0.9, na.rm = T),
                   l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                                l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                                l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                                votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                                l1.topinher	= quantile(macro$l1.topinher, 0.5, na.rm = T),
                                countryfnetherlands = 1
)												

equalLTDScenHtax = data.frame(l1.tradeallyadv = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T),  
                              l1rural.ineq  = quantile(macro$l1rural.ineq, 0.1, na.rm = T),
               #    l1.lta.ineq = quantile(macro$l1.tradeallyadv, 0.1, na.rm = T)*quantile(macro$l1rural.ineq, 0.1, na.rm = T),
                              l1.e2 = quantile(macro$l1.e2, 0.5, na.rm = T),
                              l1.trade = quantile(macro$l1.trade, 0.5, na.rm = T),
                              l1.rgdppcTenKs = quantile(macro$l1.rgdppcTenKs, 0.5, na.rm = T),
                              votetax2 = quantile(macro$votetax2, 0.5, na.rm = T),
                              l1.topinher	= quantile(macro$l1.topinher, 0.5, na.rm = T),
                              countryfnetherlands = 1
)		

scenHtax4corners = list(equalLTAScenHtax, unequalLTAScenHtax, equalLTDScenHtax, unequalLTDScenHtax)

Sim4htax <- dynsim(obj = htax.preferred.asldv, ldv = 'l1.topinher',
                  scen = scenHtax4corners, n = 20)                
dynsimGG(Sim4htax, ylab = "Top inheritance tax rate", color = 'Set1', leg.labels = c("Equal, high LTA", "Unequal, high LTA", "Equal, low LTA", "Unequal, low LTA"))

dev.print(device = pdf, file = "figureA4d.pdf", width = 7.5, height = 4)   
dev.off()

###### Table A6: Taxes as a share of GDP
size.empty = lm(Central_t_y ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                + as.factor(country_name) + as.factor(year), data = macro)

size.preferred = lm(Central_t_y ~ l1.tradeallyadv + l1rural.ineq + l1.lta.ineq
                    + l1.e2 + l1.trade + l1.rgdppc + votetax2
                    + as.factor(country_name) + as.factor(year), data = macro)


stargazer(size.empty, size.preferred
          , omit = c("year", "country", "Constant")
          , type = "text", style = "ajps", omit.stat = c("ser", "f")
          , add.lines = list(c("Country fixed effects", rep("Y", 2)), c("Year fixed effects", rep("Y", 2)))
          , digits = 2
          , title = "The determinants of domestic indirect tax revenues as a share of total tax revenue, 1871-1913"
          , order = c("tradeally", "rural.ineq", "lta.ineq", "rgdppc", "votetax", "e2", "trade")
          , covariate.labels = c("Labour trade advantage (LTA)", "Inequality", "LTA:Inequality", "GDP per capita", "Vote-tax link", "Economic franchise", "Trade")
          , dep.var.labels = "Taxes/GDP"
)


### 3. The British Case in Comparative Context
## Figure A5

ineq19001910 = c(by(macro$l1rural.ineq[which(macro$year >= 1900 & macro$year <= 1910)], macro$Country_name[which(macro$year >= 1900 & macro$year <= 1910)], mean, na.rm = T))[c(2:4, 6:13)]
ineq18701889 = c(by(macro$l1rural.ineq[which(macro$year < 1900 )], macro$Country_name[which( macro$year < 1900)], mean, na.rm = T))[c(2:4, 6:13)]

shadesOfGrey <- colorRampPalette(c("grey0", "grey100"))
threeGreys <- shadesOfGrey(3)
candp = cbind(ineq18701889, ineq19001910)
par(mar = c(4, 8, 1, 1)+0.2)
cbp = barplot(t(candp[order(candp[,2]),]), horiz = T, las = 2, beside= T, col = c(threeGreys[3], threeGreys[2]), xlab = "Average inequality", xlim = c(0, 0.7))
legend(0.5, cbp[6], fill = c(threeGreys[3], threeGreys[2]), legend = c("1870-1899", "1900-1910"), bty = "n", cex = 0.8)

dev.print(device = pdf, file = "figureA5.pdf",  width = 5.5, height = 6.5)
dev.off()

### 4. Coding the Historical Record
## Table A8: industrialisation by county/borough
load("trade_redist_brit_constit.RData") ## data is called britcon

tti1 = table(britcon$d.type, britcon$industrial.f)
tti2 = prop.table(tti1, 1)

boroughs = c(tti1[1,], sum(tti1[1,]))
share.boroughs = c(round(tti2[1,], digits = 2), NA)
counties = c(tti1[2,], sum(tti1[2,]))
share.counties = c(round(tti2[2,], digits = 2), NA)

total = apply(tti1, 2, sum)
total = c(total, sum(total))

a8 = as.data.frame(rbind(boroughs, share.boroughs, counties, share.counties, total))
a8

### 5. Ease of Labour entry in 1906
## Table A10

summary(britcon$labelease)

## Figure A7
ur = unique(britcon$region2)
labunop = rep(NA, length(ur))
for(i in 1:length(ur)){
  labunop[i] = length(which(britcon$region2 == ur[i] & britcon$labelease == 2))
  }
par(mar = c(5, 8, 2, 2)+0.2)
barplot(labunop[order(labunop)], names.arg = ur[order(labunop)], horiz = 2, las = 2, col = "darkgrey")

dev.print(device = pdf, file= "figureA7.pdf",  width = 7.5, height = 3.5)
dev.off()

## Table A11: Bivariate model -- Labour candidates unopposed

britcon$lab.unopp = as.numeric(britcon$labelease == 2)

a11 = glm(lab.unopp ~ trade.gen.f + double + d.type + lb.socec + industrial.f + contest1900 + cvs1900.2+ cvs1900.22 + scotland,
          data = britcon,
        family = binomial(link = "logit"))

a11.t.vcov<-cluster.vcov(a11, britcon$ctype, force_posdef = TRUE)
a11.coef = coeftest(a11, vcov = a11.t.vcov)

stargazer(a11.coef
          , style = "ajps"
          , title="Bivariate model: Where labor candidates run unopposed."
          , label = "tab:uk_ease_ft_unopp"
          , order = c("trade.gen.f1", "trade.gen.f0", "double", "typec", "socecB", "socecC", "industrial.f1", "industrial.f2","contest", "cvs", "scotland")
          , covariate.labels = c("Free trade", "Neutral trade", "Double member", "County", "Mixed class", "Working class", "Industrial", "Part-industrial","Election 1900: Liberal", "Election 1900: Conservative", "Election 1900: Neither", "Cons. vote 1900", "Cons. vote 1900 squared","Scotland")
          , notes.append = T
          , dep.var.labels = "Dependent variable: unopposed labor run."
          , notes = c("Standard errors clustered by constituency type", "Omitted categories are: protectionist, upper class socioeconomic,", "contested by Liberals and Conservatives in 1900, non-industrial", "single-member, borough.")
         , type = "text"
          , digits = 2
)

## 6. Votes on the People's Budget

## Table A13: Descriptive statistics for Pelling and the People's Budget

summary(britcon$lb.socec)
summary(britcon$trade.constit.f)
summary(britcon$trade.gen.f)
summary(britcon$industrial.f)
summary(britcon$d.type)
summary(as.factor(britcon$military))
summary(as.factor(britcon$aye))

## Figure A8: Votes on the Third Reading of the 1909 Finance Bill

## vote numbers
labour.n = c(length(which(britcon$aye == 1 & britcon$party2 == "Labour")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Labour")), length(which(britcon$aye == 0 & britcon$party2 == "Labour")))
liberal.n = c(length(which(britcon$aye == 1 & britcon$party2 == "Liberal")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Liberal")), length(which(britcon$aye == 0 & britcon$party2 == "Liberal")))
liberal.u.n = c(length(which(britcon$aye == 1 & britcon$party2 == "Liberal Unionist")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Liberal Unionist")), length(which(britcon$aye == 0 & britcon$party2 == "Liberal Unionist")))
unionist.n= c(length(which(britcon$aye == 1 & britcon$party == "Unionist")), length(which(is.na(britcon$aye) == T & britcon$party == "Unionist")), length(which(britcon$aye == 0 & britcon$party == "Unionist")))
conservative.n= c(length(which(britcon$aye == 1 & britcon$party == "Conservative")), length(which(is.na(britcon$aye) == T & britcon$party == "Conservative")), length(which(britcon$aye == 0 & britcon$party == "Conservative")))
irish.n = c(length(which(britcon$aye == 1 & britcon$party2 == "Irish Nationalist Party")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Irish Nationalist Party")), length(which(britcon$aye == 0 & britcon$party2 == "Irish Nationalist Party"))) +
  c(length(which(britcon$aye == 1 & britcon$party2 == "Irish Parliamentary Party")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Irish Parliamentary Party")), length(which(britcon$aye == 0 & britcon$party2 == "Irish Parliamentary Party"))) +
  c(length(which(britcon$aye == 1 & britcon$party2 == "Irish Unionist")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Irish Unionist")), length(which(britcon$aye == 0 & britcon$party2 == "Irish Unionist")))+
  c(length(which(britcon$aye == 1 & britcon$party2 == "Nationalist Party (Ireland)")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Nationalist Party (Ireland)")), length(which(britcon$aye == 0 & britcon$party2 == "Nationalist Party (Ireland)")))+
  c(length(which(britcon$aye == 1 & britcon$party2 == "Independent Nationalist")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Independent Nationalist")), length(which(britcon$aye == 0 & britcon$party2 == "Independent Nationalist")))+
  c(length(which(britcon$aye == 1 & britcon$party2 == "")), length(which(is.na(britcon$aye) == T & britcon$party2 == "")), length(which(britcon$aye == 0 & britcon$party2 == "")))

other.n = c(length(which(britcon$aye == 1 & britcon$party2 == "Independent Liberal")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Independent Liberal")), length(which(britcon$aye == 0 & britcon$party2 == "Independent Liberal"))) +
  c(length(which(britcon$aye == 1 & britcon$party2 == "Independent Unionist")), length(which(is.na(britcon$aye) == T & britcon$party2 == "Independent Unionist")), length(which(britcon$aye == 0 & britcon$party2 == "Independent Unionist")))

## vote percentages
labour = labour.n/sum(labour.n)
liberal = liberal.n/sum(liberal.n)
liberal.u = liberal.u.n/sum(liberal.u.n)
unionist = unionist.n/sum(unionist.n)
conservative = conservative.n/sum(conservative.n)
irish = irish.n/sum(irish.n)
other = other.n/sum(other.n)

vote = rbind(liberal, labour, liberal.u, conservative, unionist)
rownames(vote) = c("Liberal", "Labour", "Liberal \n Unionist", "Conservative", "Unionist")
vote = vote*100
fourGreys <- shadesOfGrey(4)

par(mar = c(4,6,4,4)+0.3)

plot(seq(0,100, length = 100), seq(0,6,length=100), col = "transparent", xaxt = "n", yaxt = "n", xlab = "", ylab = "")
mtext(rownames(vote), side = 2, at = c(1,2,3,4,5), las = 1, line = 1)

segments(x0 = rep(0, 5), x1 = vote[,1], y0 =c(1,2,3,4,5), y1 = c(1,2,3,4,5), lwd = 4, lend = 1)

segments(x0 = rep(100, 5), x1 = 100-vote[,3], y0 =c(1,2,3,4,5), y1 = c(1,2,3,4,5), lwd = 4, lend = 1, col = fourGreys[3])
axis(1)
mtext("Share of MPs voting 'Aye'", adj = 0, side = 1, line = 2.5)

axis(3, at = c(0, 20, 40, 60, 80, 100), labels = c(100, 80, 60, 40, 20, 0), col.axis = fourGreys[3])
mtext("Share of MPs voting 'No'", adj = 1, side = 3, line = 2.5, col = fourGreys[3])

mtext(side = 4, at = c(1,2,3,4,5), c(sum(liberal.n), sum(labour.n), sum(liberal.u.n), sum(conservative.n), sum(unionist.n)), line = 1, las = 2)
mtext(side = 4, at = 6, text = "No. of \n MPs", line = 1, las=2)

dev.print(device = pdf, file = "figureA8.pdf", width = 5.5, height = 3.5)
dev.off()

## Figure A9 -- Expected Values, Party Model

m2f.type = glm(aye ~ trade.gen.f + military + industrial.f + d.type + lb.socec + party2, family=binomial(link = "logit"), data = britcon, x=T)
m2f.t.vcov<-cluster.vcov(m2f.type, britcon$ctype)
m2f.t.coef = coeftest(m2f.type, m2f.t.vcov)

m2f.sim <- rmvnorm(1000, m2f.type$coef, m2f.t.vcov)

# Set the values of the independent variables
freetrade.2libx <- c(1, #intercept
                     0, # not neutral
                     1, # free trade
                     0, # not military
                     0, # not part-industrial
                     0, # not industrial
                     0, # not a county 
                     0, # not mixed socioecon
                     0, # not lower socioecon 
                     0, # not Independent Liberal
                     0, # not Labour
                     1, # Liberal
                     0 # not Liberal Unionist
                     ) 

protection.2libx <- c(1, #intercept
                         0, # not neutral
                         0, # not free trade --> omitted category is protectionist
                         0, # not military
                         0, # not part-industrial
                         0, # not industrial
                         0, # not a county 
                         0, # not mixed socioecon
                         0, # not lower socioecon 
                         0, # not Independent Liberal
                         0, # not Labour
                         1, # Liberal
                         0 # not Liberal Unionist
) 

freetrade.2conx <- c(1, #intercept
                     0, # not neutral
                     1, # free trade
                     0, # not military
                     0, # not part-industrial
                     0, # not industrial
                     0, # not a county 
                     0, # not mixed socioecon
                     0, # not lower socioecon 
                     0, # not Independent Liberal
                     0, # not Labour
                     0, # not Liberal
                     0 # not Liberal Unionist --> omitted party is Conservative/Unionist
                    ) 

protection.2conx <- c(1, #intercept
                      0, # not neutral
                      0, # not free trade --> omitted category is protectionist
                      0, # not military
                      0, # not part-industrial
                      0, # not industrial
                      0, # not a county 
                      0, # not mixed socioecon
                      0, # not lower socioecon 
                      0, # not Independent Liberal
                      0, # not Labour
                      0, # not Liberal
                      0 # not Liberal Unionist --> omitted party is Conservative/Unionist
                      
) 

pp.2freetrade <- numeric(1000)
for(j in 1:1000) {
  pp.2freetrade[j] <- t(as.matrix(freetrade.2libx)) %*% m2f.sim[j,]
}

pe.2ft <- mean(inv.logit(pp.2freetrade))
lo.2ft <- quantile(inv.logit(pp.2freetrade), .025)
hi.2ft <- quantile(inv.logit(pp.2freetrade), .975)

pp.2protection <- numeric(1000)
for(j in 1:1000) {
  pp.2protection[j] <- t(as.matrix(protection.2libx)) %*% m2f.sim[j,]
}

pe.2pr <- mean(inv.logit(pp.2protection))
lo.2pr <- quantile(inv.logit(pp.2protection), .025)
hi.2pr <- quantile(inv.logit(pp.2protection), .975)


pp.2con.freetrade <- numeric(1000)
for(j in 1:1000) {
  pp.2con.freetrade[j] <- t(as.matrix(freetrade.2conx)) %*% m2f.sim[j,]
}

pe.2conft <- mean(inv.logit(pp.2con.freetrade))
lo.2conft <- quantile(inv.logit(pp.2con.freetrade), .025)
hi.2conft <- quantile(inv.logit(pp.2con.freetrade), .975)

pp.2con.protection <- numeric(1000)
for(j in 1:1000) {
  pp.2con.protection[j] <- t(as.matrix(protection.2conx)) %*% m2f.sim[j,]
}

pe.2conpr <- mean(inv.logit(pp.2con.protection))
lo.2conpr <- quantile(inv.logit(pp.2con.protection), .025)
hi.2conpr <- quantile(inv.logit(pp.2con.protection), .975)

par(mar = c(4,6,1,1)+0.2)
bpep = barplot(height = c(pe.2ft, pe.2pr, pe.2conft, pe.2conpr), horiz = T, xlim = c(-0.05, 1.05), xlab = "Expected probability of `aye'", col = c("#a6cee3", "#a6cee3","#b2df8a", "#b2df8a"), names.arg = c("Free trade \n constituency", "Protectionist  \n constituency", "Free trade  \n constituency", "Protectionist \n constituency"), las = 1 )
segments(lo.2ft, bpep[1], hi.2ft, bpep[1], lty=3)
segments(lo.2pr, bpep[2], hi.2pr, bpep[2], lty=3)

points(x = c(lo.2ft,  hi.2ft), y = c(bpep[1],bpep[1]), pch = 18)
points(x = c(lo.2pr,  hi.2pr), y = c(bpep[2],bpep[2]), pch = 18)

segments(lo.2conft, bpep[3], hi.2conft, bpep[3], lty=3)
segments(lo.2conpr, bpep[4], hi.2conpr, bpep[4], lty=3)

points(x = c(lo.2conft,  hi.2conft), y = c(bpep[3],bpep[3]), pch = 18)
points(x = c(lo.2conpr,  hi.2conpr), y = c(bpep[4],bpep[4]), pch = 18)

legend(0.4, bpep[4], fill = c("#b2df8a", "#a6cee3"), legend = c("Conservative/ Unionist", "Liberal"), bty = "n", yjust = 0.5)

dev.print(device = pdf, file = "figureA9.pdf", width = 6.5, height = 5.5)
dev.off()

